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APPENDIX A 
COUNTERACTING GEOMETRIC 
DISTORTIONS IN WATERMARKING 



# include 
# include 
/include 
#include 
/include 



<stdio.h> 

<stdlib.h> 

<math.h> 

<gl/gl.h> 

<gl/device.h> 
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FREOS.H 



#define DC_VALUE_RATIO 0.2 

5 long number_f reqs = 4; 

double freq[10] = { 
1.9, 
1.654, 

10 1.476, 
1.227 

}; 



double phase [10] = { 
15B 2.111, 
^0 0.0765, 
1.32, 
1.0, 
2.38, 

m °- 73 < 

in 3-o, 
= i.i 

}; 
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STAMP_IT.C 

/* this program explores the specifics of stamping an image with 'digital 
reticles' for the purpose of 
5 determining registration between a given image and the digimarc signatures, 

i.e., the registration 

requirements for finding the scale and rotation elements of the ubiquitous 
snowy image patterns */ 

10 



/include "pinecone.h" 
/include "freqs.h" 

is 

^ /define DISPLAY_IT 1 

/define PI 3,141592653589 
jp /define SQRT2_2 0.707106781186 

#define INITIAL_SCALE 10.0 

m 

til unsigned char *img; 
s _ long xdim,Ydim; 

Colorindex *disp; 
;u Colorindex *scale_disp; 
?g Colorindex *pdisp; 

: 

jj^ int helplines =2; 
char *help[ ] = 
{ 

30 "usage: stamp_it filename xdim ydim channels \n", 

w \n stamp_it modifies the input image to include digital reticles and 
outputs filename. stamped " 

}; 

35 

int load_scale_disp{void) { 
long i f j; 

pdisp = scale_disp+20; 
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for( j=0;j<15; j++){ 

for(i*l;i<20;i++){ 

*pdisp = (unsigned char) 255; 
pdisp +* 20; 

} 

pdisp 20; 

} 

pdisp = scale_disp+20+ (35*400) ; 
for(j=0;j<15; j++){ 

for(i=l;i<20;i++){ 

*pdisp = (unsigned char)255; 

pdisp += 20; 

} 

pdisp += 20; 

} 

return ( 0 ) ; 
} 

int draw_scale( double value) { 
long i; 
double dtemp; 

value *=10.0; 

\ 

memcpy (disp, scale_disp, 40000) ; 
dtemp = 20,0 * value; 
if(dtemp > 399«0)dtemp=399.0; 
pdisp * &disp[ 15*400 + ( int ) dtemp ] ; 
for(i=0;i<20;i++) { 

*pdisp - (unsigned char) 255; 

pdisp += 400; 

> 

rectwr ite ( 0 , 0 , 399 , 49 , disp) ; 

return ( 0 ) ; 
} 
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double calculate_change< double distance , double scale , long which ){ 
long i; 

double value=DC_VALUE_RATIO * scale; 

for { i=0 ; i<number_f reqs ; i++ ) { 

value +« scale * sin ( freq[i] * distance + phase[i] + 
(double) which * PI ); 
} 

if (value < 0.0)value = 0.0; 



return (value) ; 
} 



main( argc, argv ) 
int argc ; 
char *argv[] ; 

{ 

long i , j , go=l , button , channels , which , count ; 
long temp, increment , last_middle; 

unsigned char tmp, *pimg, *pimgl, *img_out, *img_r, *img_b; 
char string[80], outfile[80]; 
FILE *inf; 

double change, current_scale = INITIALJSCALE; 
long gid_img,gid_stamped,gid_scale; 



if(argci=5) { 

for( j=0; j<helplines; j++)fprintf { stderr, n %s n , help(jj) 
exit( 1 ) ; 

} 



xdim = atoi(argv[2] ) ; 

ydim = atoi(argv[3 ] ) ; 

channels = atoi(argv[4) ) ; 

if (channels i= 1 && channels 1=3) { 

fprintf( stderr, tt stamp_it : channels must equal 1 for B/W 
3 for color\n H ) ; 

exit( 1 ) ; 



APPENDIX A, Page 93 



WYCtdte 43578.p« 9/25/95 at 3) 



-94- 



} 

if (DISPLAY_IT) disp - calloc (xdim*ydim, sizeof (Color index) ) ; 

else disp = calloc (40000, sizeof (Colorindex) ) ; 

scale_disp » calloc( 40000, sizeof (Colorindex) ) ; 

img = calloc (xdim*ydim, sizeof (unsigned char) ) ; 

img_out = calloc (xdim*ydim, sizeof (unsigned char) ) ; 

if{ tdisp !| !scale_disp \\ i img || I img_out ){ 

fprintf( stderr, H stamp_it : can not allocate space\n M ) ; 

exit( 1 ) ; 

} 

if (channels ~ 3){ 

img_r = calloc (xdim*ydim, sizeof (unsigned char) ) ; 
img_b = calloc (xdim*ydim, sizeof (unsigned char) ) ; 
if ( I img_r | | ! img_b ) { 

fprintf( stderr, "stamp_it : can not allocate space\ 

exit( 1 ) ; 

} 

} 

/* read in binary data into array */ 
inf = fopen(argv[l] , "r") ; 
if(iinf) { 

fprintf (stderr, "stamp_it: can't open %s\n w , argv[ 1] ) ; 
exit (1) ; 

} 

if (channels == 3){ 

f read ( img_r, sizeof (unsigned char) ,xdim*ydim, inf) ; 
f read( img, sizeof (unsigned char) ,xdim*ydim, inf) ; 
f read (img_b, sizeof (unsigned char) ,xdim*ydim, inf) ; 

} 

else f read ( img, sizeof (unsigned char ) ,xdim*ydim, inf ) ; 
fclose(inf ) ; 

/* flip it */ 
pimg = img ; 
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10 



pimgl * &img[xdim* (ydim-1) ] ; 
for(i=0;i<(ydim/2) { 

for { j=0; j<xdim; { 

tmp » *pimg; 

*(pimg++) = *pimgl; 

*(pimgl++) - tmp; 

} 

pimgl -= (2*xdim); 

} 



if (DISPLAY_IT) { 

foreground ( ) ; 
pref size(xdim,ydim) ; 
IS* gid_img = winopen ( "Original" ) ; 

?Q gflush{); 

?y pdisp - disp; 

pimg = img; 

^02 for(i=0;i<(ydim*xdim) ;i++) { 

fri *{pdisp++) = ( Color index )* (pimg++ ) ; 

> 

\ 3 rectwr ite (0,0, xdim-1 , ydim-1 , disp ) ; 

2j5f pref size(xdim,ydim) ; 

m gid_stamped - winopen { "Stamped Image..."); 

i T rectwr ite (0,0, xdim-1 , ydim-1 , disp) ; 

load_scale_disp( ) ; 
30 pref size(400, 50) ; 



gid_scale ~ winopen ( "Rough scaling..."); 



} 



35 

/* main visual feedback loop */ 
lastjmiddle = 0; 
while (go) { 
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/* first diagonal */ 
for(i=0;i<(xdim+ydim-l) ;i++) { 
which - 0; 

/* calculate addition or subtraction to image */ 
5 change » calculate_change( (double) i * SQRT2_2 , 

current_scale, which ) ; 

if ( i < xdim ) { 

pimg = &img{ i*xdim] ; 
pimgl = &img_out [ i*xdim] ; 
10 count = i+1; 

} 

else { 

pimg = &img[xdim*ydim - ydim - xdim + i + 1]; 
pimgl = &img_out [xdim*ydim - xdim - ydim + i + 

count = xdim + ydim - i - 1; 

*B > 

HI for{ j=0; j<count; j++) { 

temp = (long)*pimg + ( long) { change+0. 5) ; 
2<£ if (temp > 255) temp = 255; 

*pimgl ~ (unsigned char) temp; 
g Pimg ~= (xdim-1); 

\j pimgl -= (xdim-1); 

m } 
ifj } 

V* /* second diagonal */ 

:T for (i=0;i<{xdim+ydim-l) ;i++) { 

which = 1 ; 

/* calculate addition or subtraction to image */ 
30 change - calculate_change ( (double) (xdim - 1 - i) * 

SQRT2_2 , current_scale, which ) ; 

if ( i < xdim ) { 

pimgl = &img_out [xdim - i - 1]; 
count = i+1; 

35 } 

else { 

pimgl - &img_out[ { i-xdim+1 ) *xdim] ; 
count = xdirn+ydim-i-1; 

} 
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for( j=0; j<count; j++) { 

temp = {long)*pimgl + ( long) (change+O . 5) ; 
if (temp > 255} temp = 255; 
*pimgl = (unsigned char) temp; 
pimgl += (xdim+1); 

} 

} 



if (DISPLA¥_IT) { 

pdisp = disp; 

pimg = img_out; 

for ( i=0; i<(ydim*xdim) ; i++) { 

*(pdisp++) = { Co lor index )* (pimg++) ; 

} 

winset (gid_stamped) ; 

rectwrite(0, 0,xdim-l r ydim-i, disp) ; 

pdisp = disp; 
pimg ~ img; 

for (i=0;i<(ydim*xdim) ;i++) { 

*{pdisp++) = (Coiorindex)*(pimg++) ; 

} 

winset ( gid_img ) ; 

rectwrite(0,0,xdim-l, ydim-1, disp) ; 



winset (gid_scale) ; 
draw_scale (current_scale) ; 



button=0; 
while ( I button) { 

if( getbutton(LEFTMOUSE) ){ 

current_scale *= 1.15; 

button = 1; 

last_middle = 0; 

} 

if( getbutton(RIGHTMOUSE) ){ 

current scale *= 0.85; 
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button = 1; 
iast_middle =0; 

} 

if( getbutton(MIDDLEMOUSE) ){ 

while( getbutton(MIDDLEMOUSE) ); 
if( last_middle ){ 
button = 1; 
go=0; 

} 

last_middle = 1; 

} 

} 

} 



/* now re-sign output image because of slight changes to master key 

/* flip output image */ 
pimg = img_out; 

pimgl = &img_out [xdim* (ydim-1 ) ] ; 
for (i=0;i<(ydim/2) ;i++) { 

for ( j =0 ; j<xdim; j ++ ) { 

tmp =s *pimg; 

*(pimg++) = *pimgl; 

*(pimgl++) = tmp; 

} 

pimgl -~ (2*xdim); 

} 

/* write out signed image */ 
sprintf (outf ile, **%s. stamped" , argv{ 1 ] ) ; 
inf = f open (outf ile, "w") ; 
if(iinf) { 

fprintf ( stderr, "stamp_it : can't open %s\n" , outf ile) ; 
exit(l) ; 

} 

if (channels == 3){ 

fwrite( img_r, sizeof (unsigned char) ,xdim*ydim, inf) ; 
fwrite( img_out, sizeof (unsigned char) , xdim*ydim, inf) ; 
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fwrite(img_b,sizeof {unsigned char) ,xdim*ydim, inf ) ; 

} 

else fwrite( img_out, sizeof (unsigned char ), xdim*ydim, inf ) ; 
f close ( inf ) ; 

/* free and clean up */ 

if (DISPLAY_IT) gexit(); 

free<img) ; 

free ( img_out ) ; 

free(disp) ; 

f ree{ scale_disp) ; 

if( channels = 3){ 

f ree(img_r ) ; 

free { img_b ) ; 

} 
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REGISTR.C 

/* this program is a companion to stamp_it, wherein it is given a suspect 
image and it attempts to determine where 

the cross-hatch pattern resides within the suspect image , thereby determining 
the scale, roation and offset 
of the suspect image */ 

#include "pinecone. h" 
/include "freqs.h" 

/define DISPLAY_IT 1 
/define DISP_POW 0.1 

/define MOV_AV 21 /* keep this odd please */ 

#define MAGGJTHRESHOLD 1.7 

/define ANGLE_I NCREMENT (PI/360.0) 

/define START_SCALE_ZONE 50 

/define SCALERS TART 0.5 

/define SCALE_STOP 2.0 

/define SCALE_STEP 0.001 

/define MAX_BLOCK_SIZE 5 

/define PI 3.141592653589 

unsigned char *img; 

long xdim,ydim; 

Colorindex *disp; 

Colorindex *pdisp; 

long bits, f f tdim, f ft_size; 

float *ar,*ai; 

float wr[ 10000] ,wi[ 10000] ; 

int helplines = 2; 
char *help(] = 
{ 

"usage: register filename xdim ydim channels \n w , 
"\n register looks at the input image for digital reticles and 
displays registered result H 

}; 
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int shift_array{ float *array,int dim){ 
int i r j; 

int dim2 ~ dim/2 f 
int offset = dim2*dim + dim2; 
5 float *pl,*p2, ftmp; 



for ( i-0 ; i<dim2 ; i++ ) { 

pi = &array[i*dim] ; 
p2 = fcarray [of f set+i*dim] ; 
10 f or ( j =0 ; j <dim2 ; j ++ ) { 

ftmp = *pl; 
*pl = *p2; 
*p2 « ftmp; 
pl++;p2++; 

sp offset ~ dim2*dim; 

I U for ( i=0 ; i<d im2 ; i++ ) { 

■^H pi = &array(dim2+i*dim} ; 

2Qf p2 = fcarray [of f set+i*dim] ; 

jj^ for{j=0;j<dim2; j++){ 

ftmp = *pl; 

p *pl = *p2; 

iQ *p2 = ftmp; 

ijjj pl++;p2++; 

[| } 

li } 

return (0) ; 

} 

30 

int block_f i Iter ( float *array,int dim) { 
int i,j,k,l,i2; 
float buffer [2 J [4096] ; 



35 for(i=0;i<dim;i++)buf fer[0] [i] = array(i]; 

for{i=l;i<(dim-l) ;i++) { 
i2 = i%2; 

buffer[i2] [0] = array [ i*dim] ; 

buf fer[ 12] [dim-1] = array [ i*dim+dim-l] ; 
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(j+i)]; 



for(j-l;j<(dim-l); j++){ 

buffer[i2][ j]=0.0; 
for{k=-l;k<2;k++) { 

for (1=-1;1<2 ;!++){ 

buf fer[i2] [ j]+= array [( i+k) *dim + 

} 

} 

buffer{i2] [jJ/=9.0; 

} 

i2 = (i-l)%2; 

for ( j=0; j<dim; j++) array [ (i-l)*dim + j] = buf f er [ i2 ] [ j ] ; 

} 

12 = (i-l)%2; 

for{ j=0; j<dim; j++)array[ (i-l)*dim + j] = buf f er [ i2 ] [ j ] ; 
return ( 0 ) ; 



/* assumes fftdim arrays ar and ai * / 
double get_mag (double x, double y) { 

long xoff,yoff; 

float *par,*pai; 

double xdist , ydist , cf , revalue , i_value , value=0 . 0 ; 

xof f ~ (long)x; 

yoff ~ (long)y; 

xdist = x - (double)xof f ; 

ydist - y - { double )yof f ; 

par = &ar[yoff * fftdim + xoff]; 
pai = &ai(yoff * fftdim + xoff]; 

cf = (1.0 - xdist) * (1.0-ydist); 
revalue = (double) * <par++) ; 
i_value = (double) * (pai++) ; 
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value * cf * sqrt ( r_value* revalue + i_value*i_value) ; 
cf = xdist * (1.0-ydist); 
revalue * (double) *par; 
i_value = (double) *pai; 
5 value +* cf * sqrt (r_value* revalue + i_value*i_value) ; 

par +* (fftdim-1); 
pai +» (fftdim-1); 

cf * (1.0 - xdist) * ydist; 
10 r_value » (double) * {par ++) ; 

i_value * (double) * (pai++) ; 

value +=s cf * sqrt (revalue* revalue + i_value*i_value) ; 
cf = xdist * ydist; 
revalue ~ (double) *par; 
l&s i_value = (double) *pai; 

value +» cf * sqrt (r_value*r_value + i_value* i_value) ; 

iM' return (value) ; 

P? } 
2C| 



#5 

] J main ( argc , argv ) 
1^ int argc ; 

char *argv[ ] ; 

{ 

30 long i,j,k,go,channels, count, center; 

long dim, angle_int , total_angles , top_candidate ; 

unsigned char tmp , *pimg , *pimgl , * img_out , * img_r , * img_b ; 

char string [80], outfile{80); 

FILE *inf; 
35 long gid_img; 

double 

x,y,coss, sinn, angle, mag,magg[ 5000] ,magg_ma( 5000] , angle_vector [ 10000] ; 

double radius , dtmp, dscale , grey_dif f , highest , f rac , highest_scale, scale ; 
float *par,*pai; 
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if{argcl=5) { 

for( j=0; j<helplines; j++)fprintf ( stderr, "%s M , help[j]) ; 
exit< 1 ) ; 

} 

xdim = atoi(argv[2] ) ; 

ydim = atoi ( argv[ 3 ] ) ; 

channels = atoi(argv[4] ) ; 

if (channels i= 1 && channels 1=3) { 

fprintf ( stderr, "register : channels must equal 1 for B/W or 
3 for color\n" ) ; 

exit{ 1 ) ; 

} 

/* find the next power of two equal to or higher than the highest 
input dimension */ 

if (xdim > ydim)dim = xdim; 
else dim = ydim; 
fftdim - 1; go = 1; bits = 0; 
while ( go ) { 

if ( dim > fftdim ) { 
fftdim*=2; 
bits++; 

} 

else go = 0; 

} 

if (bits > 12) { 

fprintf ( stderr, "recognize : sorry, this particular program 
only accepts 4K images and less\n H ) ; 
exit( 1 ) ; 

} 

fft_size = fftdim * fftdim; 

disp = calloc( f f t_size, sizeof (Color index) ) ; 
img =s calloc(xdira*ydim, sizeof (unsigned char) ) ; 
ar = calloc(f ft_size, sizeof (float) ) ; 
ai = calloc(f ft_size, sizeof (float ) ) ; 
if( Idisp j| iimg [j iar || »ai ){ 

fprintf ( stderr, "register : can not allocate spaceXn" ) ; 
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exit( 1 ) ; 

} 

if (channels == 3){ 

5 img_r = calloc(xdim*ydim, sizeof (unsigned char) ) ; 

img_b - calloc(xdim*ydim, sizeof (unsigned char) ) ; 
if( iimg_r \\ 1 img Jd ){ 

fprintf( stderr, "register : can not allocate space\n" 

) ; 

10 exit( 1 ) ; 

} 

} 

/* read in binary data into array */ 
|& inf = fopenfargvfll/'r"); 

\Q if (i inf) { 

sQi fprintf (stderr, "register : can't open %s\n" , argv[ 1 ) ) ; 

jU exit(l); 

:;5 > 

?S if (channels « 3){ 

y4 f read (img_r, sizeof (unsigned char) , xdim*ydim, inf ) ; 

f read ( img, sizeof (unsigned char ), xdim*ydim, inf ) ; 
Q f read ( img_b, sizeof (unsigned char ), xdim*ydim, inf ) ; 

53 > 

23 else fread( img, sizeof (unsigned char) ,xdim*ydim, inf ) ; 

fclose(inf); 

/* flip it */ 
pimg = img; 
30 pimgl = &img[xdim* (ydim-1 ) ] ; 

for(i=0;i<(ydim/2) ;i++) { 

for ( j=0; j<xdim; j++) { 
tmp = *pimg; 
*(pirag++) = *pimgl; 
35 *(pimgl++) - tmp; 

} 

pimgl -= (2*xdim); 

} 
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/* copy image buffer into ar */ 
par = ar; 
pimg - img; 
for(i=0;i<ydim;i++) { 

for ( j =0 ; j <xdim; j ++ ) { 

*(par++) = ( float )*(pimg++); 

} 

par += (f ftdim-xdim) ; 

} 

/* 2d fft, in place */ 

printf ( n \nforward fft... \n tt ); 

f ft2d(ar, ai,bits,0, wr,wi) ; 

printf ("done \n * ) ; 

shift_array (ar, f ftdim) ; 
shift_array(ai, f ftdim) ; 
center - f ftdim/ 2; 

/* 

block_f ilter (ar , f ftdim) ; 
block_f ilter (ai, f ftdim) ; 

*/ 

if (DISPLAY_IT) { 

foreground ( ) ; 
prefsize(f ftdim, f ftdim) ; 
gid_img - winopen( "yahh" ) ; 
gf lush( ) ; 

dtmp = 

(double) ar[center*f ftdim+center+1]* (double) ar [ center*f ftdira+center+1] 
dtmp += 

(double) ai[ center *f ftdim+center+1] * (double) ai[center*f ftdim+center+1] 
dscale = pow(dtmp, DISP_POW) ; 

pdisp = disp; 
par = ar; 
pai = ai; 

for(i=0;i<(f ftdim* f ftdim) ;i++) { 



APPENDIX A, Page 106 



WYCate 43578. p* 9/25/95 (313) 



- 107 - 



dtmp = pow( (*par * *par + *pai * *pai) , 

DISP_POW ) ; 

dtmp *= {255.0 / dscale); 
if (dtmp>255.0)dtmp = 255.0; 
5 *(pdisp++) = (Colorindex) < dtmp ); 

par++;pai++; 

} 

rectwr ite { 0 , 0 , f f tdim- 1 , f f tdim- 1 , disp ) ; 

} 

10 

/* now search for the gross rotation axes */ 

for (angle = 0. 0, angle_int = 0;angle<PI/2;angle+=ANGLE_INCR£MENT, 
angle_int++) { 

coss = cos (angle); 
If 2 sinn = sin (angle); 

^0 /* fill radial vector */ 

5^ for ( i=0 ; i< { f f tdim/2 ) ; i++ ) { 

jj - ; 

radius = (double) i; 
Tj x = (double) center - radius * coss; 

2 te y = (double)center - radius * sinn; 

Iff mag ~ get_mag (x, y ) ; 

x - (double) center + radius * sinn; 
^ y = (double)center - radius * coss; 

sil ma 9 += getjnag(x,y) ; 

maggfi] = mag; 

ii > 

\& /* create moving average */ 

magg_ma[MOV_AV/2] = 0.0; 
for(i=0;i<MOV_AV;i++) { 
30 magg_ma[M0V_AV/2] += magg[i]; 

} 

magg_ma[MOV_AV/2] /- ( { double )MOV_AV ); 
f or ( i= ( MOV_AV/ 2 ) + 1 ; i< ( f f tdim/ 2 ) - { MOV_AV/ 2 ) - 1 ; i++ ) { 
magg_ma [ i ] - magg jcna [ i-1 J ; 
35 magg_ma[i] -= ((magg[i - (MOV_AV/2) - 

1 ] ) / { double ) MOV_AV) ; 

magg_ma[i] += ( <magg[ (MOV_AV/2 ) + i] )/ (double) MOVJVV) ; 

} 
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/* within prescribed 'scale zone', calculate final number for 

this angle */ 

angle_vector[angle_int ] = 0,0; 

for(i=START_SCALE_ZONE;i<(fftdim/2) - (MOV_AV/2) -l;i++){ 
mag » magg [ i J / rnagg_ma [ i ] ; 
if( mag > MAGG_THRE S HOLD ){ 

mag -= MAGG_THRESHOLD; 

angle_vector [angle_int ] +* (mag*mag) ; 

} 

} 

} 

total_angles = angle_int; 

/* sort out the T0P_CAND IDATES and find which has the best match on 
absolute scale */ 

/* chhose the highest angle_vector number for starters */ 
highest - 0.0; 

for ( angle_int=0 ; angle_int<total_angles ; angle_int++ ) { 
if ( angle_vector [ angle_int ] >highest ) { 

highest = angle_vector [ angle_int ] ; 
top_candidate = angle_int; 

} 

} 

printf ("\n\n tilt from original found = %d (top_candidate/2 ) - 

45); 

coss ~ cos ( ANGLE_INCREMENT * ( double ) top_candidate) ; 
sinn = s in ( ANGLE_ INCREMENT * (double) top_candidate) ; 
/* fill radial vector for this angle */ 
for(i=0;i<(f ftdim/2) ;i++) { 
radius = (double) i; 

x = (double) center - radius * coss; 

y = (double) center - radius * sinn; 

mag = get_mag (x, y ) ; 

x = (double) center + radius * sinn; 

y = (double) center - radius * coss; 

mag get_mag(x,y) ; 

magg{i] = mag; 
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} 

/* create moving average */ 
maggjcna[MOV_AV/2] = 0.0; 
for { i=0 ; i<MOV_AV; i++ ) { 

magg_ma [ MO V_AV / 2 ] +« magg [ i ] ; 

} 

magg_ma(MOV_AV/2] /= ( ( double )MOV_AV ); 
for ( i= (MOV_AV/2 ) +1 ; i< ( f f tdim/2 ) -(MOV_AV/2 ) -1; i++ ) { 
magg_ma[i] = magg_ma[ ; 

maggjna[i] — ((magg[i - (MOV_AV/2) - 1 ])/ (double )MOV_AV) ; 
magg_ma(i] +* < (magg[ (MOV_AV/2) + i] )/ (double) MOV_AV) ; 

} 

/* now slide the scale and find the highest point */ 
highest = 0,0; 

for(scale = SCALERS TART ; scale < SCALE_STOP;scale+=SCALE_STEP) { 
mag = 0,0; 

f or ( j =0 ; j <number_f reqs ; j ++ ) { 

radius * scale * freq[j] * (double) fftdim / PI / 2.0 
if( ( int) radius <= (l+MOV_AV/2) jj (int ) (radius+1) > 
((fftdim/2) - <MOV_AV/2) - 1) ); 

else { 

frac = radius - (double) { (int) radius); 
mag += ( 1 . 0-f rac) * (magg[ (int) radius ] / 

magg_ma[ (int) radius ] ) ; 

mag += frac*(magg[ (int) (radius+1) ] / 

magg_ma[ ( int) (radius+1) ]); 

} 

} 

if (mag > highest) { 

highest = mag; 
highest_scale = scale; 

} 

} 

printf ("\n\nscale found = %lf \n M , 1.0/highest_scale) ; 

/* now find the exact offset and orietnation, i.e. if it is flipped, 

etc. */ 



APPENDIX A, Page 109 



WYCaik* 43578. p. OO) 



-110- 



/* display the result at correct rotation and pixel size */ 
sleep{1000) ; 

5 

/* free and clean up */ 
if (DISPLAY_IT) gexit{); 
f ree(img) ; 
free(disp) ; 
10 free{ar); 

f ree(ai) ; 

if( channels = 3){ 

free{img_r) ; 
f ree( img_b) ; 

m > 

^0 } 
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